Greenhouse experiment Stephan Grassl - microbiome analysis
Yuri Pinheiro Alves de Souza
2023-06-21
Introduction
Terrestrial Ecology Department, Technical University of Munich
In a greenhouse full factorial experiment, we investigated the effects of prothioconazole and fluoxastrobin fungicide irrigation, of inoculation of the pathogenic fungus Drechslera teres and of barley cultivars on plant growth and grain yield variables, and particularly on barley root endophytic communities. Our experimental design was an approach to real agriculture practice, where generally fungicides are irrigated prophylactic before occurrence of pathogenic fungi. Additionally, aphid abundance has been measured. 60 root samples with all treatment combinations were chosen for DNA-extraction, PCR amplification and metabarcoding. For beta diversity analysis, multivariate statistics had been performed to reduce the large data sets into few principal coordinates in a MDS-CAP ordination. We hypothesized that the fungicide treatment will affect alpha and beta diversity in bacterial and fungal endophyte communities. It was found, that the fungicide and Drechslera treatments had no significant effect on alpha and beta diversity of the endophytes, but barley cultivar showed to be a significant driver of beta diversity of bacterial and fungal barley root endophyte communities.
Bacteria
Setting envriroment
loading packages
##Loading packages
library(phyloseq)
library(ggplot2)
library(readxl)
library(reshape2)
library(dplyr)
library(microbiome)
library(DESeq2)
library(writexl)
library(microeco)
library(vegan)
library(microbiomeutilities)
library(lme4)
library(lmerTest)
library(nlme)
library(metagMisc)
library(gridExtra)
library(microeco)
library(file2meco)
library(ggalluvial)
library(ggcorrplot)
library(ggpubr)
library(car)
library(DT)
## Setting theme
theme_set(theme_classic())Loading data
setwd("C:/Users/yurip/Desktop/COMI/Colaborations/Stephan/second_round")
otu_mat<- as.data.frame(read_excel("asvs.xlsx"))
tax_mat<-
as.data.frame(read_excel("taxonomy.xlsx",
na = "NA",
col_names =
c("ASV", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus",
"Species"))) %>%
purrr::map_df(~gsub(
pattern = "metagenome|uncultured|unidentified|Unknown",
replacement = NA,
.x)) %>%
mutate_if(is.character, stringr::str_trim) %>%
mutate(Kingdom = ifelse(is.na(Kingdom),
"U. Kingdom",
Kingdom),
Phylum = coalesce(Phylum,
ifelse(grepl("^U.", Kingdom),
Kingdom,
paste("K_", Kingdom))),
Class = coalesce(Class,
ifelse(grepl("^U.", Phylum),
Phylum,
paste("P_", Phylum))),
Order = coalesce(Order,
ifelse(grepl("^U.", Class),
Class,
paste("C_", Class))),
Family = coalesce(Family,
ifelse(grepl("^U.", Order),
Order,
paste("O_", Order))),
Genus = coalesce(Genus,
ifelse(grepl("^U.", Family),
Family,
paste("F_", Family))),
Species = coalesce(Species,
ifelse(grepl("^U.", Genus),
Genus,
paste("G_", Genus)))) %>%
tibble::column_to_rownames("ASV") %>%
filter(Kingdom %in% "Bacteria" &
!grepl("Chloroplast", Genus) &
!grepl("Mitochondria", Genus) &
!Order %in% "Chloroplast")
samples_df <- as.data.frame(read_excel("metadata_stephan2.xlsx"))
otu_names=otu_mat$OTUS
row.names(otu_mat) <- otu_names
otu_mat <- otu_mat[,2:61]
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select (-sample)
samples_df$Block <- factor(samples_df$Block, levels = c("Block1", "Block2", "Block3", "Block4","Block5" ))
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
class(otu_mat) <- "numeric"
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
stephan <- phyloseq(OTU, TAX, samples)
stephan## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6275 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 6275 taxa by 7 taxonomic ranks ]sample_names(stephan)## [1] "Weisser.1_S34" "Weisser.3_S35" "Weisser.5_S36" "Weisser.7_S37"
## [5] "Weisser.9_S38" "Weisser.11_S39" "Weisser.13_S40" "Weisser.15_S41"
## [9] "Weisser.17_S42" "Weisser.19_S43" "Weisser.21_S44" "Weisser.23_S45"
## [13] "Weisser.25_S46" "Weisser.27_S47" "Weisser.29_S48" "Weisser.31_S49"
## [17] "Weisser.33_S50" "Weisser.35_S51" "Weisser.37_S52" "Weisser.39_S53"
## [21] "Weisser.41_S54" "Weisser.43_S55" "Weisser.45_S56" "Weisser.47_S57"
## [25] "Weisser.49_S58" "Weisser.51_S59" "Weisser.53_S60" "Weisser.55_S61"
## [29] "Weisser.57_S62" "Weisser.59_S63" "Weisser.61_S64" "Weisser.63_S65"
## [33] "Weisser.65_S66" "Weisser.67_S67" "Weisser.69_S68" "Weisser.71_S69"
## [37] "Weisser.73_S70" "Weisser.75_S71" "Weisser.77_S72" "Weisser.79_S73"
## [41] "Weisser.81_S74" "Weisser.83_S75" "Weisser.85_S76" "Weisser.87_S77"
## [45] "Weisser.89_S78" "Weisser.91_S79" "Weisser.93_S80" "Weisser.95_S81"
## [49] "Weisser.97_S82" "Weisser.99_S83" "Weisser.101_S84" "Weisser.103_S85"
## [53] "Weisser.105_S86" "Weisser.107_S87" "Weisser.109_S88" "Weisser.111_S89"
## [57] "Weisser.113_S90" "Weisser.115_S91" "Weisser.117_S92" "Weisser.119_S93"rank_names(stephan)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"sample_variables(stephan)## [1] "Plant_number" "Cultivar"
## [3] "Fungicide_treatment" "Drechslera_treatment"
## [5] "Aphid_treatment" "Block"
## [7] "Cultivar_and_Drechslera" "Fungicide_and_Drechslera"
## [9] "Aphid_count" "Total_Seed_Biomass"
## [11] "Dry_Weight"rank_names(stephan)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"sample_variables(stephan)## [1] "Plant_number" "Cultivar"
## [3] "Fungicide_treatment" "Drechslera_treatment"
## [5] "Aphid_treatment" "Block"
## [7] "Cultivar_and_Drechslera" "Fungicide_and_Drechslera"
## [9] "Aphid_count" "Total_Seed_Biomass"
## [11] "Dry_Weight"Rarefaction curve and normalization
mat <- t(otu_table(stephan))
class(mat) <- "matrix"
rarecurve(mat, step=50, cex=0.5, label = TRUE,
xlab = "Reads", ylab = "ASVs")
#Normalize number of reads in each sample using median sequencing depth.
total = median(sample_sums(stephan))
standf = function(x, t=total) round(t * (x / sum(x)))
stephan.rarefied = transform_sample_counts(stephan, standf)Diversity analyzes
Alpha diversity
alpha_bac=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
measures=c("Shannon")) + geom_boxplot() +
facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
scale_color_manual(values = c("#90EE90", "#228822")) +
labs(title = "Alpha diversity (Shannon) - Bacteria")
alpha_bac
alpha_bac2=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
measures=c("Simpson")) + geom_boxplot() +
facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
scale_color_manual(values = c("#90EE90", "#228822")) +
labs(title = "Alpha diversity (Simpson) - Bacteria")
alpha_bac2
## Calculation diversity metrics
diversity=estimate_richness(stephan.rarefied)
shannon=diversity$Shannon
observed=diversity$Observed
observed.log <- log(observed)
evenness=shannon/observed.log
diversity$evennes=evenness
richness=as.data.frame(diversity)
DT::datatable(richness, class = "cell-border stripe", rownames = T,
filter = "top", editable = TRUE, extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = '300px',
scrollCollapse = TRUE
))Beta diversity
### PCoA and Constrained Analysis of Principal Coordinates (CAP)
ps1V4.rel_bray <- phyloseq::distance(stephan.rarefied, method = "bray") # CAP ordinate
cap_ord <- ordinate(physeq = stephan.rarefied,
method = "CAP",
distance = ps1V4.rel_bray,
formula = ~ Aphid_count + Total_Seed_Biomass +
Dry_Weight)
scree.cap <- plot_scree(cap_ord, "Scree Plot for MCs in Constrained Analysis of Principal Coordinates - CAPSCALE")
print(scree.cap)
cap_plot <- plot_ordination(physeq = stephan.rarefied,
ordination = cap_ord,
color = "Cultivar",
shape = "Fungicide_treatment",
axes = c(1,2)) +
geom_point(aes(colour = Cultivar), size = 3) +
geom_point(size = 3) +
scale_color_manual(
values = c( "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724"))
#Now add the environmental variables as arrows
cap_ord## Call: capscale(formula = distance ~ Aphid_count + Total_Seed_Biomass +
## Dry_Weight, data = data)
##
## Inertia Proportion Rank
## Total 18.9706062 1.0000000
## Constrained 1.3895020 0.0732450 3
## Unconstrained 17.5853700 0.9269799 56
## Imaginary -0.0042658 -0.0002249 1
## Inertia is squared Bray distance
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3
## 0.6710 0.4388 0.2797
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 3.194 1.321 1.047 0.874 0.732 0.663 0.591 0.536
## (Showing 8 of 56 unconstrained eigenvalues)arrowmat <- vegan::scores(cap_ord, display = "bp")
# Add labels, make a data.frame
arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1,
yend = CAP1,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3 * CAP1,
y = 1.3 * CAP1,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cap_plot +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf,
color = "black",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf,
show.legend = FALSE
) +
stat_ellipse(type = "t") + ggtitle("CAP_Plot - Bacteria")
Taxonomic composition
Box plot
mycols <- c("yellow", "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")
pn <- plot_taxa_boxplot(stephan.rarefied,
taxonomic.level = "Family",
top.otu = 9,
group = "Cultivar_and_Drechslera",
title = "top 9 taxa at family level - Bacteria",
keep.other = FALSE,
add.violin= FALSE,
group.order = c("ACCO Drechs_NO", "ACCO Drechs_YES",
"KLAR Drechs_NO", "KLAR Drechs_YES" ,
"PLAN Drechs_YES", "PLAN Drechs_NO"),
group.colors = mycols
)
box_plot_bacteria = pn + theme_biome_utils() + theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
theme(legend.position = "bottom")
box_plot_bacteria
### Exporting data most abundant taxa abundances for linear models ###
top_nine_bacteria <- subset_taxa(stephan.rarefied, Family %in% c("Cellvibrionaceae" , "Comamonadaceae",
"Enterobacteriaceae", "Oxalobacteraceae" ,
"Pseudomonadaceae", "Rhizobiaceae",
"Rhodanobacteraceae", "Xanthobacteraceae",
"Xanthomonadaceae"))
top_nine_bacteria## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1298 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 1298 taxa by 7 taxonomic ranks ]top_nine_melted=psmelt(top_nine_bacteria)
write_xlsx(top_nine_melted,"top_nine_bacteria.xlsx")Heatmap
meco_dataset <- phyloseq2meco(stephan.rarefied)
t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 10, groupmean = "Fungicide_and_Drechslera")
t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 20)
t1$plot_heatmap(facet = "Cultivar", xtext_keep = FALSE, withmargin = FALSE)
Statistical analyzes
Permanova
stephan_bray <- phyloseq::distance(stephan, method = "bray")
stephandf <- data.frame(sample_data(stephan))
adonis2(stephan_bray ~ Fungicide_treatment*Drechslera_treatment*Cultivar, data = stephandf)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = stephan_bray ~ Fungicide_treatment * Drechslera_treatment * Cultivar, data = stephandf)
## Df SumOfSqs R2 F
## Fungicide_treatment 1 0.3135 0.01637 0.9617
## Drechslera_treatment 1 0.2633 0.01375 0.8077
## Cultivar 2 1.0527 0.05498 1.6146
## Fungicide_treatment:Drechslera_treatment 1 0.2451 0.01280 0.7519
## Fungicide_treatment:Cultivar 2 0.4920 0.02570 0.7546
## Drechslera_treatment:Cultivar 2 0.5335 0.02786 0.8183
## Fungicide_treatment:Drechslera_treatment:Cultivar 2 0.5986 0.03126 0.9181
## Residual 48 15.6482 0.81727
## Total 59 19.1470 1.00000
## Pr(>F)
## Fungicide_treatment 0.464
## Drechslera_treatment 0.754
## Cultivar 0.009 **
## Fungicide_treatment:Drechslera_treatment 0.858
## Fungicide_treatment:Cultivar 0.929
## Drechslera_treatment:Cultivar 0.830
## Fungicide_treatment:Drechslera_treatment:Cultivar 0.624
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Linear mixed effect models
richness$Block = samples_df$Block
richness$Drechslera_treatment = samples_df$Drechslera_treatment
richness$Fungicide_treatment = samples_df$Fungicide_treatment
richness$Sample = samples_df$Plant_number
richness$Cultivar = samples_df$Cultivar
table(richness$Fungicide_treatment,richness$Drechslera_treatment)##
## Drechs_NO Drechs_YES
## Fungicide_NO 15 15
## Fungicide_YES 15 15m1<-lme(Shannon~Drechslera_treatment*Fungicide_treatment*Cultivar, ,
random=~1|Block,
na.action=na.exclude,
data=richness)
anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 44 1158.7360 <.0001
## Drechslera_treatment 1 44 1.6947 0.1998
## Fungicide_treatment 1 44 0.0191 0.8908
## Cultivar 2 44 0.8794 0.4222
## Drechslera_treatment:Fungicide_treatment 1 44 0.5659 0.4559
## Drechslera_treatment:Cultivar 2 44 0.9959 0.3776
## Fungicide_treatment:Cultivar 2 44 0.8806 0.4217
## Drechslera_treatment:Fungicide_treatment:Cultivar 2 44 0.0467 0.9544plot(m1, main = "Residuals Shannon alpha diversity")
m2<-lme(Simpson~Drechslera_treatment*Fungicide_treatment*Cultivar, ,
random=~1|Block,
na.action=na.exclude,
data=richness)
anova(m2)## numDF denDF F-value p-value
## (Intercept) 1 44 10405.612 <.0001
## Drechslera_treatment 1 44 1.177 0.2838
## Fungicide_treatment 1 44 0.129 0.7212
## Cultivar 2 44 0.631 0.5366
## Drechslera_treatment:Fungicide_treatment 1 44 0.070 0.7921
## Drechslera_treatment:Cultivar 2 44 2.003 0.1471
## Fungicide_treatment:Cultivar 2 44 1.206 0.3090
## Drechslera_treatment:Fungicide_treatment:Cultivar 2 44 1.526 0.2286plot(m2, main = "Residuals Simpson alpha diversity")
Fungi
Loading data
setwd("C:/Users/yurip/Desktop/COMI/Colaborations/Stephan/ITS")
otu_mat<- as.data.frame(read_excel("asvs.ITS.xlsx"))
tax_mat<-
as.data.frame(read_excel("taxonomy.ITS.xlsx",
na = "NA",
col_names =
c("ASV", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus",
"Species"))) %>%
purrr::map_df(~gsub(
pattern = "metagenome|uncultured|unidentified|Unknown",
replacement = NA,
.x)) %>%
mutate_if(is.character, stringr::str_trim) %>%
mutate(Kingdom = ifelse(is.na(Kingdom),
"U. Kingdom",
Kingdom),
Phylum = coalesce(Phylum,
ifelse(grepl("^U.", Kingdom),
Kingdom,
paste("K_", Kingdom))),
Class = coalesce(Class,
ifelse(grepl("^U.", Phylum),
Phylum,
paste("P_", Phylum))),
Order = coalesce(Order,
ifelse(grepl("^U.", Class),
Class,
paste("C_", Class))),
Family = coalesce(Family,
ifelse(grepl("^U.", Order),
Order,
paste("O_", Order))),
Genus = coalesce(Genus,
ifelse(grepl("^U.", Family),
Family,
paste("F_", Family))),
Species = coalesce(Species,
ifelse(grepl("^U.", Genus),
Genus,
paste("G_", Genus)))) %>%
tibble::column_to_rownames("ASV")
samples_df <- as.data.frame(read_excel("metadata_stephan.xlsx"))
otu_names=otu_mat$ASV
row.names(otu_mat) <- otu_names
otu_mat <- otu_mat[,2:61]
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select (-sample)
samples_df$Block <- factor(samples_df$Block, levels = c("Block1", "Block2", "Block3", "Block4","Block5" ))
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
class(otu_mat) <- "numeric"
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
stephan <- phyloseq(OTU, TAX, samples)
stephan## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1082 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 1082 taxa by 7 taxonomic ranks ]rank_names(stephan)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"sample_variables(stephan)## [1] "Plant_number" "Cultivar"
## [3] "Fungicide_treatment" "Drechslera_treatment"
## [5] "Aphid_treatment" "Block"
## [7] "Cultivar_and_Drechslera" "Fungicide_and_Drechslera"
## [9] "Aphid_count" "Total_Seed_Biomass"
## [11] "Dry_Weight"Rarefaction curve and normalization
mat <- t(otu_table(stephan))
class(mat) <- "matrix"
rarecurve(mat, step=50, cex=0.5, label = TRUE,
xlab = "Reads", ylab = "ASVs")
#Normalize number of reads in each sample using median sequencing depth.
stephan.2 <- subset_samples(stephan, !(Plant_number %in% c("27")))
total = median(sample_sums(stephan.2))
standf = function(x, t=total) round(t * (x / sum(x)))
stephan.rarefied = transform_sample_counts(stephan.2, standf)Diversity analyzes
Alpha diversity
stephan.rarefied <- subset_samples(stephan.rarefied, !(Plant_number %in% c("27"))) ## I still need to put the new sample 27 in the dataset
alpha_fung=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
measures=c("Shannon")) + geom_boxplot() +
facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
scale_color_manual(values = c("#90EE90", "#228822")) +
labs(title = "Alpha diversity (Shannon) - Fungi")
alpha_fung
alpha_fung2=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
measures=c("Simpson")) + geom_boxplot() +
facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
scale_color_manual(values = c("#90EE90", "#228822")) +
labs(title = "Alpha diversity (Simpson) - Fungi")
alpha_fung2
## Calculation diversity metrics
diversity=estimate_richness(stephan.rarefied)
shannon=diversity$Shannon
observed=diversity$Observed
observed.log <- log(observed)
evenness=shannon/observed.log
diversity$evennes=evenness
richness=as.data.frame(diversity)
DT::datatable(richness, class = "cell-border stripe", rownames = T,
filter = "top", editable = TRUE, extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = '300px',
scrollCollapse = TRUE
))Beta diversity
### PCoA and Constrained Analysis of Principal Coordinates (CAP)
ps1V4.rel_bray <- phyloseq::distance(stephan.rarefied, method = "bray") # CAP ordinate
cap_ord <- ordinate(physeq = stephan.rarefied,
method = "CAP",
distance = ps1V4.rel_bray,
formula = ~ Aphid_count + Total_Seed_Biomass +
Dry_Weight)
scree.cap <- plot_scree(cap_ord, "Scree Plot for MCs in Constrained Analysis of Principal Coordinates - Fungi")
print(scree.cap)
cap_plot <- plot_ordination(physeq = stephan.rarefied,
ordination = cap_ord,
color = "Cultivar",
shape = "Fungicide_treatment",
axes = c(1,2)) +
geom_point(aes(colour = Cultivar), size = 3) +
geom_point(size = 3) +
scale_color_manual(
values = c( "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724"))
#Now add the environmental variables as arrows
cap_ord## Call: capscale(formula = distance ~ Aphid_count + Total_Seed_Biomass +
## Dry_Weight, data = data)
##
## Inertia Proportion Rank
## Total 11.8754 1.0000
## Constrained 1.3749 0.1158 3
## Unconstrained 12.0069 1.0111 34
## Imaginary -1.5064 -0.1269 24
## Inertia is squared Bray distance
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3
## 0.9305 0.3480 0.0964
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 5.857 1.536 1.080 0.539 0.478 0.401 0.373 0.249
## (Showing 8 of 34 unconstrained eigenvalues)arrowmat <- vegan::scores(cap_ord, display = "bp")
# Add labels, make a data.frame
arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1,
yend = CAP1,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3 * CAP1,
y = 1.3 * CAP1,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cap_plot +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf,
color = "black",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf,
show.legend = FALSE
) +
stat_ellipse(type = "t") + ggtitle("CAP_Plot - Fungi")
Taxonomic composition
Box plot
mycols <- c("yellow", "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")
tax_table(stephan.rarefied)[,colnames(tax_table(stephan.rarefied))] <- gsub(tax_table(stephan.rarefied)[,colnames(tax_table(stephan.rarefied))],pattern="O_ C_ P_ K_ k__Fungi",replacement="Unidentified fungi")
pn <- plot_taxa_boxplot(stephan.rarefied,
taxonomic.level = "Family",
top.otu = 9,
group = "Cultivar_and_Drechslera",
title = "top 9 taxa at family level - Fungi",
keep.other = FALSE,
add.violin= FALSE,
group.order = c("ACCO Drechs_NO", "ACCO Drechs_YES",
"KLAR Drechs_NO", "KLAR Drechs_YES" ,
"PLAN Drechs_YES", "PLAN Drechs_NO"),
group.colors = mycols
)
box_plot_fungi = pn + theme_biome_utils() + theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
theme(legend.position = "bottom")
box_plot_fungi
top_nine_fungi <- subset_taxa(stephan.rarefied, Family %in% c("Chaetomiaceae" , "Helotiaceae",
"Naviculisporaceae", "Nectriaceae" ,
"Podosporaceae", "Unidentified fungi", "Pleurotheciaceae", "O_C_P_p__Ascomycota", "O_C_P_p__Chytridiomycota"))
top_nine_fungi## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 738 taxa and 59 samples ]
## sample_data() Sample Data: [ 59 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 738 taxa by 7 taxonomic ranks ]top_nine_melted_fungi=psmelt(top_nine_fungi)
write_xlsx(top_nine_melted_fungi,"top_nine_fungi.xlsx")Heatmap
meco_dataset <- phyloseq2meco(stephan.rarefied)
t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 10, groupmean = "Fungicide_and_Drechslera")
t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 20)
t1$plot_heatmap(facet = "Cultivar", xtext_keep = FALSE, withmargin = FALSE)
Statistical analyzes
Permanova
stephan_bray <- phyloseq::distance(stephan, method = "bray")
stephandf <- data.frame(sample_data(stephan))
adonis2(stephan_bray ~ Fungicide_treatment*Drechslera_treatment*Cultivar, data = stephandf)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = stephan_bray ~ Fungicide_treatment * Drechslera_treatment * Cultivar, data = stephandf)
## Df SumOfSqs R2 F
## Fungicide_treatment 1 0.1234 0.00982 1.0979
## Drechslera_treatment 1 0.1537 0.01223 1.3679
## Cultivar 2 6.1188 0.48684 27.2221
## Fungicide_treatment:Drechslera_treatment 1 0.2027 0.01613 1.8035
## Fungicide_treatment:Cultivar 2 0.2036 0.01620 0.9057
## Drechslera_treatment:Cultivar 2 0.2158 0.01717 0.9600
## Fungicide_treatment:Drechslera_treatment:Cultivar 2 0.1559 0.01241 0.6937
## Residual 48 5.3946 0.42921
## Total 59 12.5685 1.00000
## Pr(>F)
## Fungicide_treatment 0.305
## Drechslera_treatment 0.199
## Cultivar 0.001 ***
## Fungicide_treatment:Drechslera_treatment 0.126
## Fungicide_treatment:Cultivar 0.481
## Drechslera_treatment:Cultivar 0.436
## Fungicide_treatment:Drechslera_treatment:Cultivar 0.682
## Residual
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Linear mixed effect models
samples_df = samples_df %>% filter(Plant_number != 27)
richness$Block = samples_df$Block
richness$Drechslera_treatment = samples_df$Drechslera_treatment
richness$Fungicide_treatment = samples_df$Fungicide_treatment
richness$Sample = samples_df$Plant_number
richness$Cultivar = samples_df$Cultivar
table(richness$Fungicide_treatment,richness$Drechslera_treatment)##
## Drechs_NO Drechs_YES
## Fungicide_NO 15 14
## Fungicide_YES 15 15m1<-lme(Shannon~Drechslera_treatment*Fungicide_treatment*Cultivar, ,
random=~1|Block,
na.action=na.exclude,
data=richness)
anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 43 1755.2338 <.0001
## Drechslera_treatment 1 43 0.0000 0.9961
## Fungicide_treatment 1 43 0.0005 0.9822
## Cultivar 2 43 0.2101 0.8113
## Drechslera_treatment:Fungicide_treatment 1 43 0.1546 0.6961
## Drechslera_treatment:Cultivar 2 43 0.1601 0.8525
## Fungicide_treatment:Cultivar 2 43 1.7716 0.1823
## Drechslera_treatment:Fungicide_treatment:Cultivar 2 43 1.1499 0.3262plot(m1, main = "Residuals Shannon alpha diversity")
m2<-lme(Simpson~Drechslera_treatment*Fungicide_treatment*Cultivar, ,
random=~1|Block,
na.action=na.exclude,
data=richness)
anova(m2)## numDF denDF F-value p-value
## (Intercept) 1 43 2252.2469 <.0001
## Drechslera_treatment 1 43 0.3411 0.5622
## Fungicide_treatment 1 43 0.0576 0.8115
## Cultivar 2 43 0.2144 0.8079
## Drechslera_treatment:Fungicide_treatment 1 43 0.0617 0.8050
## Drechslera_treatment:Cultivar 2 43 0.2077 0.8132
## Fungicide_treatment:Cultivar 2 43 1.2332 0.3015
## Drechslera_treatment:Fungicide_treatment:Cultivar 2 43 0.6009 0.5529plot(m2, main = "Residuals Simpson alpha diversity")